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Abstract 

The  symbol  of  the  Hessian  for  an  aeroelastic  optimization  model  problem  is  ana¬ 
lyzed.  The  flow  is  modeled  by  the  small-disturbance  full  potential  equation  and  the 
structure  is  modeled  by  an  isotropic  (von  Karman)  plate  equation.  The  cost  function 
consists  of  both  aerodynamic  and  structural  terms.  In  the  new  analysis  the  symbol  of 
the  cost  function  Hessian  near  the  minimum  is  computed.  The  result  indicates  that 
under  some  conditions,  which  are  likely  fulfilled  in  most  applications,  the  system  is 
decoupled  for  the  non-smooth  components.  The  result  also  shows  that  the  structure 
part  in  the  Hessian  is  well-conditioned  while  the  aerodynamic  part  is  ill-conditioned. 
Applications  of  the  result  to  optimization  strategies  are  discussed. 
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1  Introduction 


Lately,  there  is  a  growing  interest  in  Multidisciplinary  Design  and  Optimization  (MDO) 
[l]-[4],  An  important  problem  in  that  field  is  the  aeroelastic  optimal  design  problem  (for 
example,  [5]-[7]).  In  this  problem  there  are  two  coupled  disciplines:  aerodynamics  and 
structural  analysis.  The  problem  is  to  compute  the  aerodynamic  shape  and  structural  rigidity 
such  that  some  given  cost  function  is  minimized. 

The  purpose  of  this  paper  is  to  demonstrate  new  analysis  of  Hessians  for  MDO  problems 
on  the  above  aeroelastic  optimization  problem  and  to  draw  some  practical  conclusions.  The 
approach  is  to  consider  a  simple  model  problem  and  compute  the  symbol  of  the  Hessian  near 
the  minimum  for  the  non-smooth  frequencies.  The  Hessian  contains  curvature  information 
which  is  essential  for  the  solution  of  ill-conditioned  optimization  problems.  Hessian  symbols 
were  previously  computed  for  smoothing  predictions  in  the  development  of  multigrid  one-shot 
methods  [8]- [11]  and  lately  for  the  analysis  of  inviscid  aerodynamic  optimization  problems 
[12].  The  analysis  in  this  paper  indicates  that  for  the  non-smooth  components  the  system 
is  decoupled  under  certain  conditions,  which  are  likely  fulfilled  in  most  applications.  The 
analysis  also  shows  that  the  structures  part  in  the  Hessian  is  well- conditioned  while  the 
aerodynamics  part  is  ill-conditioned. 

One  consequence  of  this  result  is  that  if  the  decoupling  condition  holds  the  solution  of 
such  problems  can  be  achieved  in  two  stages.  In  the  first  stage,  the  MDO  approach  should 
be  taken  on  a  coarse  model;  that  is,  the  flow  and  the  structure  equations  are  considered 
simultaneously  during  the  minimization,  which  is  a  more  complex  problem  than  optimizing 
the  decoupled  individual  disciplines  problems.  In  the  second  stage,  a  refined  CFD  code  for 
the  flow  and  a  detailed  finite  element  code  for  structure  should  be  used  in  a  serial  algorithm  in 
which  the  shape  is  optimized  relative  to  aerodynamic  considerations,  followed  by  structural 
optimization  limited  to  a  given  shape  that  under  flow  conditions  it  will  twist  and  bend  to  the 
aerodynamic  optimal  shape  [13,  14].  This  approach  should  result  in  a  good  approximation 
of  the  multidisciplinary  optimal  solution. 

A  second  consequence  is  that  if  the  decoupling  condition  does  not  hold,  and  thus  the  two 
disciplines  are  coupled,  a  multidisciplinary  algorithm  is  necessary  and  a  preconditioner  for 
the  coupled  system  is  necessary  to  obtain  effective  convergence. 

The  paper  outline  is  as  follows.  In  Section  2  the  optimization  problem  is  formulated. 
In  Section  3  the  necessary  conditions  for  a  minimum  are  derived  with  the  adjoint  method 
and  their  relation  with  the  Hessian  is  discussed.  In  Section  4  the  symbol  of  the  Hessian  for 
the  non-smooth  frequencies  is  derived  by  using  local  mode  analysis.  In  Section  5  decoupling 
conditions  and  multidisciplinary  preconditioners  are  derived.  In  Section  6  applications  of 
the  result  are  finally  discussed. 


2  Problem  Formulation 

In  this  section  the  aeroelastic  analysis  problem  and  the  optimal  design  problem  are  presented. 
The  aeroelastic  analysis  problem  couples  the  full  potential  flow  equation  with  the  isotropic 
von  Karman  plate  equation  to  give  the  pressure  distribution  over  the  plate,  p,  and  the 
plate  deformation,  W,  for  a  given  plate  shape,  a,  and  rigidity  distribution,  D.  The  design 
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problem  is  to  compute  the  “best”  shape  and  structural  rigidity  so  that  a  given  cost  function 
is  minimized. 

The  cost  function  is  composed  of  aerodynamic  and  structure  parts.  The  aerodynamic 
cost  function  estimates  performance  by  measuring  the  difference,  in  L2  norm,  of  the  pressure 
distribution  from  a  desired  one.  The  structure  cost  function  gives  a  measure  of  the  structural 
weight  and  penalizes  structural  deformation. 

Since  our  interest  is  in  a  local  mode  analysis  of  the  Hessian  near  the  minimum,  we  consider 
the  small  disturbance  equations  of  flow  over  a  flat  plate. 


2.1  The  Flow  Model 

We  choose  the  full  potential  equation  as  a  model  for  the  flow.  It  approximates  inviscid 
flow  characteristics  and  is  used  in  applications  for  aerodynamic  optimal  design  (for  example, 
[15]).  For  the  analysis  of  the  cost  function’s  Hessian  in  the  vicinity  of  the  minimum  it  is 
enough  to  consider  small  perturbations  of  the  shape  from  the  optimal  solution.  The  resulting 
changes  in  the  potential  satisfy  the  steady  state  small  disturbance  full  potential  equation. 
The  geometry  is  taken  to  be  half-space  ft  =  {x,y,z  >  0),  where  the  x  axis  is  the  stream- 
wise  coordinate,  y  is  the  coordinate  perpendicular  to  the  stream  and  parallel  to  the  plate 
(span wise  direction),  and  z  is  in  the  normal  direction  to  the  plate. 

The  Aerodynamic  State  Equation: 


L<f>  =  0  z>0  ^  i\ 

B</>  =  C(a  +  W)  z  =  0  1  ; 

with  the  following  definitions  of  the  interior  operator,  L,  and  boundary  operators,  B  and  C, 

L  =  (1  —  M^dxx  +  dyy  +  dzz 

B  =  dz  (2.2) 

C  =  Uoodx 

where  Uoo  is  the  free  stream  velocity,  M  is  the  Mach  number,  with  the  following  far-field 
boundary  conditions: 

Inflow  Boundary  Condition 
Subsonic:  <f>x  =  Uoo 

Supersonic:  <j>x  =  Uoo  and  <j)  =  0  (we  assume  that  the  normal  free  stream  velocity,  Ko,  is 
zero). 

Outflow  Boundary  Condition 
Subsonic:  <f>x  =  Uoo 
Supersonic:  No  Boundary  Condition. 

The  missing  low-order  terms  in  the  boundary  condition  of  (2.1)  vanish  if  the  analysis  is 
performed  around  a  flat  shape. 


2 


2.2  The  Structural  Model 


The  structural  model  consists  of  the  isotropic  von  Karman  plate  equation  for  the  displace¬ 
ment  W  [16,  17]: 

G{D,W)  =  pC<f>  z  =  0  (2.3) 

with  the  following  definition  of  the  the  operator  G: 

G  (-D,  W)  =  ( DWXX)XX  +  (DWyy)yy  +  V  ( DWyy)  xx  +  (DWXX)yy  +  2(l  —  V )  (  D WX  y )  Xy  ,  (2.4) 

where  D  is  the  plate  rigidity  distribution,  p  is  the  flow  density  and  v  is  the  Poisson  ratio.  In 
two  space  dimensions  Eq.(2.3)  reduces  to  the  beam  equation: 

{DWXX)XX  =  pUoo(j)x  z  =  0.  (2.5) 

The  boundary  conditions  for  the  plate  are  that  it  is  clamped  at  y  =  0  (or  at  x  =  0  for  a 
beam)  and  free  at  the  other  boundaries: 

W'(®,0)  =  W'y(*,0)  =  0. 

However,  Eq.(2.3)  is  elliptic,  so  the  eifect  of  a  high-frequency  error  component  in  the  deflec¬ 
tion  W  is  local,  and  therefore  the  plate  boundary  conditions  do  not  play  a  role  in  the  local 
mode  analysis. 

2.3  The  Cost  Function  Model 

The  definition  of  the  cost  function  is  not  unique  and  depends  on  the  specific  application 
under  consideration.  In  general  the  requirement  of  the  aeroelastic  optimal  design  is  that  it 
have  maximum  aerodynamic  performance  and  minimum  structural  weight  and  deformation. 
Some  of  the  desired  features  of  the  final  design  are  in  many  cases  modeled  by  a  set  of 
inequality  constraints,  as  is  the  case  for  the  minimum  deformation  requirement.  However, 
for  the  purpose  of  this  paper  we  will  avoid  inequality  constraints  by  adding  a  term  to  the 
cost  function  which  penalizes  deformation.  In  the  following  the  different  terms  composing 
the  cost  function  are  discussed. 

•  The  Aerodynamic  Performance  Term 

A  common  aerodynamic  cost  function  is  drag  (or  drag  over  lift).  However,  in  inviscid 
aerodynamic  optimization  models  a  commonly  used  cost  function  is  pressure  matching  (for 
example,  [18]-[22]).  In  potential  models  the  pressure,  p,  is  related  to  the  potential,  <f>,  by  the 
Bernoulli  relation 

P  —  pUoo4*X- 

This  results  in  the  following  cost  function  term 

Faero  =  Jr(<j>x-  r? da 

where  da  is  an  integration  element  on  the  shape  T.  The  target  distribution,  f*  (x,y)  €  T2(r), 
is  related  to  the  desired  pressure  distribution,  p*(x,y),  by  the  relation 


•  The  Structural  Weight  Term 

Another  important  factor  in  aeroelastic  design  is  the  resulting  weight  of  the  structure.  In 
practice  the  weight  is  measured  by  the  sum  of  the  weights  of  all  the  components  composing 
the  structure.  In  plate  models  the  weight  is  related  with  the  plate  rigidity,  D ,  and  is  given  in 
the  following  table,  where  E  is  the  Young  modulus  of  elasticity,  b  and  h  are  the  cross  section 


D 

weight 

beam 

EbP  ' 
12 

Ppfr  bhdx 

plate 

Et 3 

12(1— i/2) 

pp  fr  tdxdy 

components  of  the  beam,  pP  is  the  structural  density,  and  t  is  the  plate’s  thickness. 

In  both  cases  the  weight  of  the  structure  is  proportional  to  D*  where  d  is  the  space 
dimension:  . 

pwdght  ^  j  D*da. 

•  The  Structural  Deformation  Term 

As  a  result  of  the  pressure,  p ,  exerted  on  the  plate  by  the  flow,  the  structure  will  deform 
its  shape  by  W  (bend  and  twist).  In  practice  the  structure  is  designed  so  that  the  amount  of 
deformation  will  be  constrained  not  to  exceed  some  given  limits.  In  this  model  we  account 
for  this  requirement  by  penalizing  the  deformation  which  is  measured  by  the  work  of  the 
aerodynamic  pressure  on  the  plate,  pW.  This  will  add  to  the  cost  function  the  following 
term  (using  the  Bernoulli  relation): 

pdeform  =  pUoo  j^Wda. 

2.4  The  Optimization  Problem 

We  define  the  cost  function,  F  =  F(4>,  W,  D),  to  be 

F(<f>,  W,  D )  =  7i  J(4X  -  f*fda  + 12  jf  D?da  +  73  <S>xWda  (2.6) 

where  71,72,73  are  parameters.  The  cost  function  is  a  map  from  a  function  space  to  1R. 

The  minimization  problem  is  to  find  a  shape  function,  a,  and  rigidity  distribution,  D, 
such  that  the  cost  function  is  minimized  subject  to  Eqs.(2.1)  and  (2.3).  We  assume  the 
existence  of  a  solution  for  both  the  state  equations  and  for  the  optimization  problem  (a 
rigorous  treatment  of  existence  and  uniqueness  of  solutions  is  beyond  the  scope  of  this 
paper). 

3  Adjoint  Formulation  and  the  Hessian 

In  this  section  the  necessary  conditions  for  a  minimum  are  derived  with  the  adjoint  method 
[19]-[24].  The  necessary  conditions  are  given  as  a  set  of  state  equations  (the  analysis  prob¬ 
lem),  costate  equations  (the  adjoint  problem)  and  design  equations  (optimality  conditions). 
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Then  the  relation  between  the  design  equation  residuals  and  the  Hessian  of  the  cost  function 
is  discussed.  This  relation  will  be  used  in  the  next  section  to  derive  the  Hessian’s  symbol. 


3.1  The  Necessary  Conditions  for  a  Minimum 

The  Lagrangian  is  a  functional  defined  by 

C(t,W,a,D,(,K<l)  =  F(t, W, D )  +  Jr f (B«> -  C(a  +  W))ia+ 

fa  XL<j)dQ,  +  /r  77  (g (D,  W)-pC<j>)d<T  [  ’ 

where  £  =  £(x,y),  r\  =  ri(x,y)  and  A  =  X(x,y,z )  are  Lagrange  multipliers.  The  first  order 
necessary  conditions  for  a  minimum  are  derived  by  the  requirement  that  the  first  order 
variation  of  the  Lagrangian  vanish  (this  is  known  as  the  adjoint  method  and  the  resulting 
conditions  are  known  as  the  Kuhn- Tucker  conditions). 

When  considering  the  variation  of  the  structure  state  equation  a  linearization  is  per¬ 
formed, 

G(D*  +  D,W*  +  W)  =  G(ZT,  W*)  +  Gd(W*)D  +  Gw(D*)W  +  h.o.t.  (3.2) 

where  D  and  W  are  small  perturbations  of  the  displacement  and  rigidity  from  the  optimal 
solution  W*  and  D*  respectively,  and  where  the  linearized  operators  Gj)  and  G\y  are  defined 
as  follows 

GB(W*)D  =  bxxw;x  +  Dyyw;y  +  v[bxxw;y  +  byyw;x]  +  2(1  -  v)bxyw;y 

Gw(D*)W  =  G(D*,W).  1  ’ 

Formally,  W*  and  D*  serve  as  non-constant  coefficients  in  the  linearized  structure  operator. 
In  the  following  the  costate  and  design  equations  are  given. 

Costate  Equations 

LA  =  0  z  >  0 

BA  -1-  pC 77  =  F4  z  =  0  (3.4) 

Gw(D*)r]  +  CA  =  —Fw  z  =  0 

Inflow  Boundary  Condition 
Subsonic:  A^  =  0 

Supersonic:  No  Boundary  Condition 

Outflow  Boundary  Condition 
Subsonic:  Xx  —  0 
Supersonic:  A  =  0  and  Xx  =  0 


Design  Equations 


CA  =  0  z  =  0 
Gd  +  Fd  =  0  2  =  0 


(3.5) 
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where 

^  =  -27l(<^ -/*)*- 73^4 

Fw  -  73 <t>x  (3.6) 

Fd  =  JD 

and  where  the  operators  in  the  adjoint  and  design  equations  (3.4-3.5)  satisfy 

L  =  L 

C  =  -C 
Gw  (D*)  =  Gvv(^*) 

GD(W*)  =  Gj>(W*). 

The  adjoint  boundary  operator  B  corresponds  to  the  normal  derivative,  dz,  applied  to  a 
solution  of  the  interior  costate  PDE,  A,  when  using  the  adjoint  far-held  boundary  conditions. 
We  assume  the  existence  of  a  solution  to  the  costate  equations. 


3.2  The  Relation  of  the  Hessian  with  the  Necessary  Conditions 

If  the  state  and  costate  equations  are  satisfied,  in  the  strong  sense,  then  the  variation  of  the 
Lagrangian  (3.1)  is  equal  to  the  variation  in  the  cost  function  and  is  given  by 


SC  =  f  aC \da  +  ^  D  (GD( W*)j]  +  FD)  da 


where  a  and  D  are  variations  in  the  design  variables.  Therefore  the  quantities  multiplying 
a  and  D  in  (3.7)  are  the  sensitivity  gradients  of  the  cost  function  with  respect  to  the  design 
variables,  when  computed  on  the  constraint  manifold: 

J\  =  CA  /o  o\ 

<?2  =  Gd(W*)77  +  Fd- 

The  state  and  costate  equations,  (2.1),  (2.3)  and  (3.4),  give  an  implicit  relation  between  the 
costate  variables  and  the  design  variables: 

A  =  A(o:,  D )  /o  qx 

V  =  V(*,D). 

Using  equations  (3.8)  and  (3.9)  we  can  write  the  following  relation  which  holds  near  the 
minimum  (or*  and  D*  denote  the  optimal  value  of  the  design  variables  a  and  D,  respectively): 


gi(A(o:*  +  a,  D*  +  /)))  =  H\\Oc  +  Fl^iF)  -f  h.o.t. 
g2(D*  +  D, T}(a*  +  a,D*  +  D ))  =  H21a  +  H22D  4-  h.o.t. 


(3.10) 


where  at  the  minimum 


gi(\(a*,D*))  =  g2(D*Ma*,D*))  =  0. 


We  conclude  that  on  the  constraint  manifold,  near  the  minimum,  the  Hessian  of  the  cost 
function  relates  the  errors  in  the  design  variables  with  the  residuals  of  the  design  equations 
(sensitivity  gradients).  In  the  next  section  we  will  use  this  fact  to  calculate  the  symbol  of 
the  Hessian. 
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4  Derivation  of  the  Hessian’s  Symbol 

In  the  following  section  we  compute  the  symbol  of  the  Hessian  with  local  mode  analysis. 
Hessian  symbols  were  previously  computed  for  smoothing  prediction  in  the  development 
of  multigrid  one-shot  method  [8]-[ll]  and  lately  for  the  analysis  of  inviscid  aerodynamic 
optimization  problems  [12].  In  the  following  the  local  mode  analysis  is  outlined. 

The  analysis  is  performed  in  the  vicinity  of  the  minimum  where  the  design  variables  are 
assumed  to  have  an  error  d  and  D.  We  assume  that  the  state  and  costate  equations  are 
satisfied  and  consider  the  errors  in  the  state  and  costate  variables  (<^,  W',  A,  77)  with  respect 
to  their  value  at  the  optimal  solution.  These  errors  are  assumed  to  satisfy  homogeneous 
equations  similar  to  Eqs.(2.1,3.4,  3.5),  and  a  linearization  of  Eq.(2.3).  We  then  consider  the 
high-frequency  errors  in  the  design  variables  and  compute  an  explicit  solution  of  the  problem 
in  terms  of  exponential  functions  in  a  half-space.  Then  with  a  standard  procedure  the 
problem  in  a  half-space  is  reduced  to  the  boundary.  On  the  boundary  we  study  the  mapping 
from  the  transformed  design  variables  errors  to  the  residuals  of  the  design  equations,  (g\,g2). 
The  symbol  of  this  mapping  gives  the  eigenvalues  of  the  Hessian. 

4.1  Fourier  Representation 

We  start  with  the  Fourier  representation  of  the  solution  in  a  half-space  and  perform  local 
mode  analysis.  Consider  errors  in  the  solution  of  the  form 

d(x,  y )  =  ,  U}2)ei(UlX+w^ 

D(x,y)  =  b{co1,u2)e^x+w*y\  1  '  ’ 

As  a  result,  the  errors  in  the  state  and  costate  variables  are  assumed  to  have  the  following 
form: 

Kx,y,z)  =  ^(^i,u;2,W3)ez(a'i:r+a;2j/+“3i!) 

\(x,y,z)  =  X(L0Uu}2,^zy^lX+W2y+u,3z)  (4  2} 

W(x,y)  =  W(  co1,u2)ei^x+u,2y'> 
fj(x,y)  =  rj  (uq ,  u>2 )  e*(wi*+w2») . 

Before  computing  the  relation  between  the  state  and  costate  error  symbols,  (<^,  A,  W,  77) ,  and 
the  design  error  symbols,  (a,  D),  we  reduce  the  problem  to  the  boundary  by  eliminating  u>3 
from  the  symbols  <j)  and  A. 


4.2  Reduction  to  the  Boundary 

The  reduction  to  the  boundary  is  done  by  eliminating  u3  from  the  symbol  expressions  using 
the  interior  equations.  The  following  discussion  regarding  the  choice  of  w3  was  done  in  [12] 
and  is  repeated  here  for  completeness. 

The  term  <f>  satisfies  the  interior  equation  for  <j>: 

L^(ujuu2,u3)ei{uJlX+W2y+u^  =  0. 


(4.3) 


Assuming  a  non-trivial  solution,  ^  0,  Eq.(4.3)  results  in  an  algebraic  relation  between  u>i, 
lo2  and  u> 3: 

(1  —  M2)u\  +  +  u>1  =  0.  (4.4) 

The  choice  of  lo3  should  be  done  such  that  it  will  result  in  a  physical  solution.  We  differentiate 
between  subsonic  and  supersonic  flows. 


4.2.1  Subsonic  Flow 

In  the  subsonic  regime  (M  <  1)  the  physical  solution  is  given  by 

u>3  =  1  -  M2)  +  w|, 

which  corresponds  to  decaying  solutions: 

i(x,y,z)  =  4>{ui ,  U;2)ei(Wia;+a’23/>e-(V^(1-M2)+wl)z 
A  (x,y,z)  =  A(Wl ,  u2)e^lX+^y) 

In  that  case  the  symbols  of  the  boundary  operators,  B  and  B,  are  given  by  (recall  that  B 
and  B  are  the  normal  derivatives  applied  to  a  solution  of  the  interior  state  and  costate  PDE 
respectively) 

B  —  B  =  —  (1  -  .M2)  +  wf*  (4.5) 


4.2.2  Supersonic  Flow 

We  differentiate  between  two  supersonic  cases  which  are  determined  by  a  flow  speed  denoted 
Mc  and  given  by 

Mc  — 

The  case  1  <  M  <  Mc  results  in  the  same  symbols  for  B  and  B  as  for  the  subsonic  flow  case 
(Eq.(4.5)). 

In  the  case  Mc  <  M  both  signs  of  u>3  in  (4.4)  correspond  to  physical  solutions.  The 
positive  root  correspond  to  the  characteristic  which  propagates  into  the  shape,  <j>+,  and  the 
negative  root  correspond  to  the  characteristic  which  propagates  out  of  the  shape,  </>_,  (and 
a  similar  expression  for  A): 


(4-6) 

Since  the  inflow  information  does  not  change  as  a  result  of  a  shape  perturbation,  the 
following  holds: 

<i>+(x,y,z)  =  o.  (4-7) 

In  the  same  manner  the  outflow  characteristic  of  the  adjoint  is  not  changing  as  a  result 
of  a  shape  perturbation: 

\-(x,y,z)  =  0.  (4.8) 


Therefore 


%x,y,z)  =  ^_(w! ,  W 

A(s,y,z)  =  A+  (w! ,  ujje^^+VKli-^MW . 

We  conclude  that  for  flow  speeds  Mc  <  M  the  boundary  operator  B  is  antisymmetric, 
(with  respect  to  the  adjoint  operation),  and  the  symbols  B  and  B  are  given  by 


B  =  -isJ\ul{\-M*)+ut\ 
B  =  -  M*)  +  u f\. 


In  all  flow  conditions  the  multiplication  BB  results  in  the  same  expression: 

BB  =  |uq2(l-M2)  +  w||.  (4.10) 

By  eliminating  from  the  transformed  equations  the  state  and  costate  flow  equations  can 
be  written  on  the  surface  (uq,^)  which  corresponds  to  the  boundary  (x, y). 


4.3  Treatment  of  the  Structure  Equations 

In  this  subsection  we  give  a  short  note  concerning  the  transformation  of  the  structure  state 
and  costate  equations.  The  structure  state  and  costate  equations  contain  non-constant 
coefficients  which  should  be  frozen  prior  to  the  local  mode  analysis.  The  structure  state  and 
costate  error  equations  are  given  by  (see  Eqs. (2.3, 3.2, 3.4)) 

Gd(W*)D  +  Gw(D*)W  =  PCj>  z  =  0  , 

Gw(£>*)7?  +  CA  =  -Fw  z  =  0  [  J 

where  D,  W,  4>,  fj  and  A  denote  the  error  variables,  Fw  denotes  the  error  in  Fw,  and  the 
operators  Gd  and  Gw  are  defined  in  (3.3).  Since  Eqs. (4. 11)  have  variable  coefficients,  D* 
and  W*,  it  is  necessary  to  freeze  them  around  a  point  on  the  boundary.  This  procedure  is 
justified  as  long  as  the  errors  in  the  design  variables  are  highly  oscillatory  compared  to  W * 
and  D*.  We  denote  the  values  of  W*(x,  y)  and  D*(x,y )  at  a  point  (x0,  yo)  on  the  boundary 
by  Wq  and  Dq,  respectively. 


4.4  The  Coupled  State  and  Costate  Equations  in  Fourier  Space 

In  terms  of  their  Fourier  representation  on  the  boundary,  the  state  and  costate  error  equations 
are  given  by  the  following  matrix  form: 

(  B  -C  0  0  \  /  \  / 

-pc  gw{Dq)  00^  w  -gd{w;)d 

—F<t><j>  —F$w  B  pC  A  0 

K  Fwt  0  b  bw{D*0 )  /  V  /  V  0 
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The  various  symbols  are  given  explicitly  by 


Fu  —  271^ 

Fpw  =  —*73wi 

Fw<t>  =  H3wi 

Gw{D*)  =  D*M+LolY  +  l.o.t. 

6d(ws)  =  -a f(w;„  +  *w;m)  -  +  vw;„)  - ^(2(1  -  *)W£„) 

C  —  'iUqoLO\. 


(4.13) 


Note  that  the  terms  originating  in  the  cost  function  serve  as  a  coupling  symmetric  block 
between  the  state  and  costate  systems. 


4.5  The  Symbol  of  the  Hessian 

The  design  equations  residuals,  in  the  transformed  space,  are  given  by 

9'=|N  _  (4.14) 

h  =  GdWW  +  Fdd(DI)D 

where 

Fdd  =  12(1  ~  d)  (d*0)^ 

and  the  symbols  and  g2  are  the  symbols  of  the  sensitivity  gradients  g\  and  g2l  respectively 
(see  (3.8)). 

Inverting  the  system  (4.12)  and  substituting  A  and  r\  in  the  symbol  of  the  design  residuals 
(4.14)  results  in  a  relation  between  the  residuals  of  the  design  equations  and  the  errors  in 
the  design  variables.  In  Fourier  space, 


#11  H12 

#21  #22 


(4.15) 


where  the  matrix  #;j  is  the  symbol  of  the  Hessian,  as  discussed  in  Sec.  3.2.  #n  is  the  symbol 
of  the  aerodynamic  optimization  Hessian,  #22  of  the  structural  optimization  Hessian,  and 
#12,  #21  are  the  coupling  terms.  In  the  following,  the  terms  Hij  are  given  explicitly: 


#  —  ^FFGw{Fj^Gw^  +  2  pCFjw) 

( BGW  -  C2p)(BGw  -  C2p ) 

^  _  CGpjCF^Gw  ±  BFjwGw  +  P&hw) 

( BGW  -  C2p)(BGw  -  C2p) 

^  _  CGd(CFhGw  +  BFjwGw  ±  pC2hw) 
( BGW  -  C2p)(BGw  -  C2p) 

fj  CGUFuC  +  (B  +  B)Ftw)  - 
(. BGW  -  C2p)(BGw  -  C2p) 


#21  = 


(4.16) 


(4.17) 


(4.18) 


(4.19) 


Since  Gw  is  a  fourth  order  polynomial  in  cu1)2,  G\d,  F<t»h  B  and  B  are  of  second  order,  F<i>w 
and  C  are  first  order,  and  Fp  is  of  zero  order,  the  principal  parts  of  the  above  expressions 
(the  asymptotic  limits  of  high-frequencies),  are  given  by 


=  27i  Ul 


u>i(l  -M2)  +  u>| 


(4.20) 


C2GdFh  2 7l Ul  “t  +  vW*yy)  +  ojl(W*yy  +  vWqxx)  +  2w1u;2(1  -  u)W*xy) 

12“  2lK  BBGW  "  D&  \u2(l-M2)+u2\.(^  +  u2y 

(4.21) 

#22  *  Fdd  =  7z(1J2~  ^  {D*)  ^ .  (4.22) 

Note  that  for  simplicity  we  assumed  a  complex  representation  of  the  errors,  (4.1),  and 
obtained  a  complex  Hessian  symbol.  If  considering  a  real  representation,  i.e., 

cc(x,y)  =  a{io1,uj2)ei^x+W2^  +  aconj(LJ1,u;2)e-i^x+^ 

where  aconj  is  the  complex  conjugate  of  d,  and  a  similar  expression  for  D,  then  the  resulting 
Hessian  symbol  is  real  and  symmetric. 


4.6  Discretization  and  the  Condition  Number 

In  practice  the  problem  is  solved  numerically  and  thus  discretization  is  introduced.  Therefore 
the  analysis  should  be  performed  in  the  discrete  space  and  the  Hessian  will  depend  on  the 
specific  discretization.  For  the  “ideal”  discretization,  the  symbol  of  the  Hessian  is  equal  to 
the  differential  one  with  the  substitution 

01  02 
hi  h2 

where  (h\,h2)  are  the  mesh  sizes  in  the  (x,y)  directions,  respectively,  and  where  0i  and  02 
vary  in  the  domain  [ — 7r,  7t] . 

Note  that  “high-frequencies”  are  those  which  obey  u >■  c  for  some  constant,  c,  which  is 
determined  by  the  different  parameters  in  the  problem.  In  the  discrete  space  this  corresponds 
to  0;  >>  chi.  Since  the  constant  c  is  independent  of  the  mesh-size  h,  as  the  grid  is  refined  the 
portion  of  high-frequencies  in  the  spectrum  increases  and  therefore  the  approximation  taken 
by  the  local  mode  analysis  above  is  more  accurate  for  a  larger  part  of  the  spectrum.  This  is 
not  surprising  since  as  the  grid  is  refined  its  resolution  increases  while  the  resolution  of  the 
smooth  components  remains  unchanged. 

The  maximum  eigenvalue  of  each  of  the  disciplinary  Hessians  is  estimated  by 

■^mai  =  Fa 

Unfortunately  the  lowest  eigenvalue  cannot  be  estimated  by  the  procedure  above  since  this 
is  precisely  the  spectrum  range  in  which  the  approximation  taken  by  the  local  mode  analysis 


(aq,u;2)  =  ( 
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does  not  hold.  However,  the  lowest  eigenvalue  is  a  fixed  number,  independent  of  the  mesh- 
size,  and  therefore  the  condition  number  of  the  Hessian  is  proportional  to  Amax.  For  a 
two-dimensional  flow  over  a  beam,  (u> 2  =  0),  we  get  for  the  aerodynamic  part  of  the  Hessian 

27il%*2  1 

|1  —  M2\  h2' 

We  conclude  that  the  aerodynamic  part  of  the  Hessian  is  ill-conditioned  and  its  condition 
number  is  increasing  as  the  grid  is  refined  (see  [12]  for  further  discussion).  The  structure’s 
symbol  (4.22)  approaches  a  constant,  for  the  high-frequencies,  independent  of  the  mesh-size. 
We  therefore  conclude  that  the  structural  optimization  problem  is  well-conditioned. 


(see  Eq.(4.20)) 


5  On  the  Coupling  Between  Aerodynamic  and  Struc¬ 
tural  Design 

In  the  previous  section  we  computed  explicitly  the  Hessian’s  symbol.  The  coupling  between 
the  two  disciplines,  during  optimization,  is  determined  both  by  the  off-diagonal  terms  in 
the  symbol  of  the  Hessian  and  by  the  smoothness  of  the  design  variables.  We  say  that  the 
optimization  problem  can  be  decoupled  if  (see  Eq.(4.15)): 

\HU&\  >  \Hi*D\  (5-1) 

(or  equivalently  if  |^22^i|  ^  H1292 | ) ■  If  condition  (5.1)  holds  then  the  aerodynamic  opti¬ 
mization  problem  can  be  solved  at  the  first  stage,  independent  of  the  structural  optimization 
problem,  setting  the  structural  deformation  to  zero  ( W  =  0).  Then,  at  the  second  stage, 
the  resulting  shape  and  potential  are  given  as  inputs  to  the  structural  optimization  problem 
(see  more  details  in  the  discussion).  The  solution  of  this  two-stage  approach  will  be  a  good 
approximation  of  the  multidisciplinary  solution,  in  the  high-frequencies.  We  emphasize  that 
if  condition  (5.1)  does  not  hold  then  this  procedure  will  result  in  a  poor  approximation  of 
the  multidisciplinary  optimal  solution.  In  that  case  a  multidisciplinary  preconditioner 
should  be  constructed  to  give  a  “corrected”  disciplinary  search  directions.  This  is  done  by 
transforming  Eq.(4.15)  back  to  the  PDE  level.  Note  that  if  |jff2i<*|  ^  \&22D\  holds  and 
condition  (5.1)  does  not  hold,  the  process  of  first  computing  an  optimal  structure  and  then 
computing  an  optimal  shape  is  not  well  defined.  In  the  following  we  examine  condition  (5.1) 
explicitly  for  two-  and  three-dimensional  flows  and  illustrate  the  derivation  of  a  precondi¬ 
tioner  for  a  simple  case. 


5.1  Two  Space  Dimensions 

In  a  two  dimensional  flow  over  a  beam  the  principal  part  of  the  Hessian  is  given  by  (see 
Eqs.(4.20)-(4.22)) 


H(u i  >  1)  =  2 


*nUl  _ 

|1— M2|  DJ 


w* 

wQxx 


l lUl,  ws** 

|1— M2|  D0* 

-172  DT* 


\ 

/ 


(5.2) 
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In  practice  the  curvature  of  the  deflection  in  the  stream-wise  direction,  W£xx,  is  negligible 
with  respect  to  the  rigidity,  Dq,  as  is  the  case  for  airfoils.  Setting  W£xx  =  0  in  (5.2)  results  in 
a  diagonal  matrix  which  implies  a  complete  decoupling  between  aerodynamic  and  structure 
design. 

In  case  Wqxx  is  not  negligible,  let  us  assume  the  optimal  solution  exists  and  belongs  to 
the  following  spaces 

ae  Cs(r)  and  D  G  Ck{T) 
where  s  and  k  are  integers.  The  matrix  (5.2)  implies 


Hna  oc  u>2  s 
H12D  oc  u>~k 


(5.3) 


and  therefore  the  decoupling  condition  (5.1)  is  satisfied  if 

k  >  s  -  2  (5.4) 

i.e.,  as  long  as  the  rigidity  D  is  not  much  more  “rough”  than  the  shape  a  (D  G  Cs_2(r)  or 
rougher). 

In  case  condition  (5.4)  is  not  satisfied  a  preconditioner  should  be  constructed.  The 
disciplinary  search  directions,  a  and  D,  are  solutions  of  a  PDE  which  is  obtained  by  the 
inverse  transform  of  Eq.(4.15): 


—acaxx  —  b2  a  =  eg x  -  bg2 
—acDxx  -  b2D  =  —bgi  -  a(g2)xx 


(5.5) 


with 


a  = 


_  u*o 
|1— M2| 


271^ 


w: 


°  ~  |1-M2|  D 
C  =  D 


Using  these  directions  as  the  search  directions,  in  a  multidisciplinary  optimization  algorithm, 
should  result  in  a  much  more  effective  convergence  properties. 


5.2  Three  Space  Dimensions 

In  a  three  dimensional  configuration  we  differentiate  between  the  stream- wise  and  spanwise 
directions.  Let  us  assume  that  the  curvature  of  the  deflection  in  the  stream-wise  direction 
is  negligible,  i.e.  set  Wqxx  =  0.  As  a  result  the  coupling  term  H\2  has  the  following  form: 


TT  HV*  m  27l^  Oyy  +  U2W0yy  +  ~  ^Woxy) 

i2(  j~  D*0  \0J2(l-M*)  +  L022\-(u2  +  U2y 


(5.6) 


For  the  design  of  the  structure  in  the  spanwise  direction  only,  i.e.,  freezing  the  stream-wise 
design  as  done  in  practice  for  aircraft  wing  design,  the  off-diagonal  terms  in  the  Hessian 
vanish,  (uq  =  0),  and  therefore  the  problem  is  decoupled. 
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For  errors  in  the  stream-wise  direction  only,  (i.e.  u?2  =  0),  the  off-diagonal  terms  in  the 
Hessian  reduce  to 

A,  ,.27l  Ulvw;yy 

Hu(u2  -  0)  ~  Z).|1_Af2|  • 

By  a  similar  argument  to  the  two-dimensional  case,  if  a  €  Cs (r)  and  D  €  C,s“1( r)  or 
smoother  then  the  problem  can  be  decoupled.  If  the  structure  design  variables  belong  to  a 
much  rougher  space  than  the  shape  variables,  ( D  €  Cs~2  or  rougher),  then  a  preconditioner 
should  be  applied  similar  to  Eq.(5.5)  with 


—  \1-M2{ 

l  _  tnUL  i "Oh 

0  ~  \1-M2\  D* 


The  case  a  £  Cs  and  D  £  Cs~2  needs  a  more  careful  examination.  In  this  special  case, 
the  decoupling  condition  (5.1)  implies 


\Hn\  <  |ffn| 


I  uWQyy  I 


Let  us  assume  a  wing  like  geometry  where  a  plate  of  length  L  is  clamped  at  (y  —  0)  and  free 
at  the  other  boundaries.  Assume  that  at  the  optimal  solution  the  plate  is  bent  towards  the 
tip  ( y  =  L)  as  a  quadratic  in  y, 


W*(x,y)  =  Wtip^,  (5.8) 

where  Wtip  is  the  maximum  deflection  at  the  tip  of  the  plate.  In  that  case,  the  decoupling 
condition  (5.7)  becomes 

>24(1  —  v2)vWtip 

I  EPL1  <  ’  ' 

where  E  is  the  Young  modulus  of  elasticity,  v  is  the  Poisson  ratio  and  t  is  the  plate  thick¬ 
ness.  In  case  condition  (5.9)  is  not  satisfied  then  the  problem  can  not  be  decoupled  and  a 
preconditioner  is  required. 


6  Discussion  and  Concluding  Remarks 

The  symbol  of  the  Hessian  for  aeroelastic  optimization  model  problem  was  computed  (Eqs.(4.16- 
4.19)).  The  result  indicates  that  for  the  non-smooth  components  the  system  is  decoupled 
under  certain  conditions.  In  two  dimensions  it  is  enough  to  assume  that  the  curvature  of 
the  deflection  is  negligible,  as  is  the  case  for  airfoils,  to  obtain  decoupling  (see  Sec.5.1  for 
further  discussion  of  the  decoupling  condition).  In  three  dimensions  this  requirement  is  not 
enough  and  unless  the  rigidity  is  not  much  rougher  than  the  shape  the  system  is  coupled 
(see  Sec. 5.2). 

The  result  also  shows  that  the  aerodynamic  optimization  problem  is  ill-conditioned,  and 
therefore  second  order  information  is  essential  for  efficiently  solving  this  part  of  the  problem 
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[12],  while  the  structural  optimization  problem  is  well  conditioned.  Thus,  it  is  anticipated 
that  the  number  of  optimization  iterations  required  to  solve  the  multidisciplinary  problem 
is  determined  by  the  aerodynamic  optimization  part  of  the  problem. 

We  now  discuss  the  application  of  this  result  to  optimization  strategies  for  the  solution 
of  the  problem.  We  differentiate  between  two  basic  approaches,  the  “disciplinary”  and  the 
“multidisciplinary”.  In  the  disciplinary  approach  the  solution  of  the  problem  is  divided  so 
that  one  discipline  optimization  problem  is  solved  at  each  stage,  decoupled  from  the  other 
discipline.  In  the  multidisciplinary  approach  both  the  analysis  and  optimization  solutions 
are  performed  in  a  tightly  coupled  manner.  These  two  approaches  are  now  presented  in  more 
detail. 

•  The  Disciplinary  Approach  -  Weak  Coupling 

A  common  practical  strategy  used  to  solve  large  aeroelastic  shape  optimization  prob¬ 
lems  is  the  disciplinary  approach,  i.e.,  design  the  aerodynamic  optimal  shape  to  give 
the  best  performance  and  then  design  a  minimum  weight  structure,  restricted  to  the 
aerodynamic  shape,  that  under  flow  will  twist  and  bend  to  the  aerodynamic  optimal 
shape  [13,  14]: 

The  Disciplinary  Algorithm 


1.  The  aerodynamic  shape  optimization  problem  is  solved  setting  W  =  0, 

mina  fT(<t>x  -  f  *)2  dc 
subject  to 
L<f>  =  0  z  >  0 
B(f)  =  Ca  z  =  0. 

2.  The  structure  minimum  weight  problem  is  solved  given  a  fixed 
shape  and  pressure  (potential)  distribution,  i.e. 

mine  7  /r  D* da  -f  7'  fr  <j>xW da 
subject  to 

G{D,  W)  =  pC<j)  z  =  0 
where  7  and  7'  are  parameters. 


3.  The  final  shape,  /? ,  is  computed  such  that  under  cruise 
conditions,  (given  pressure  distribution  -  p) ,  the  shape  will 
deform  into  the  aerodynamic  optimal  shape  a. 
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•  The  Multidisciplinary  Approach  -  Tight  Coupling 

Lately  there  has  been  an  effort  to  develop  new  optimization  strategies  which  couple 
the  two  disciplines  tightly  during  the  analysis  and  optimization  computation.  This  is 
known  as  the  MDO  approach  [l]-[7].  According  to  this  approach  after  each  call  to  the 
optimizer  the  analysis  and  adjoint  equations  are  relaxed,  or  solved  exactly,  depending 
on  the  feasibility  choice  (Multi-Disciplinary  Feasibility  (MDF),  Individual  Discipline 
Feasibility  (IDF)  or  All  at  Once  (A AO),  [3]). 

The  MDO  Algorithm 


The  coupled  aerodynamic  shape  and  structure  minimum  weight  optimization 
problems  are  solved  simultaneously.  In  order  to  achieve  efficient  conve¬ 
rgence  both  disciplines  gradients  play  a  role  in  each  disciplinary 
optimizer  (see  Eq.(5.5)), 

min^c  7a  fT(<f>x  -  f  *)2  dx  +  72  fr  D^dx  +  73  fr  <j>xW dx 
subject  to 

L<j>  =  0  z  >  0 

B<f>  =  C(a+W)  z  =  0 

G(D,W)  =  pC<j>  z  =  0. 

The  aim  of  the  MDO  approach  is  to  couple  a  refined  CFD  code  with  a  detailed  finite- 
element  structural  analysis  code  to  compute  the  aeroelastic  states  prior  to  each  optimization 
iteration.  The  computational  complexity  of  the  MDO  algorithm  is  much  greater  than  that 
of  the  disciplinary  algorithm  since  at  each  multidisciplinary  iteration  both  the  aerodynamic 
and  structural  optimization  problems  have  to  be  solved  (the  MDO  approach  also  introduces 
a  technical  difficulty  of  joining  together  two  large  codes).  The  result  of  this  paper  shows  that 
under  certain  conditions,  given  in  Sec. 5,  which  are  most  likely  met  for  aircraft  wing  design, 
the  system  is  decoupled  for  the  non-smooth  frequencies.  Therefore,  under  these  conditions 
the  MDO  approach  applied  on  a  fine  scale  model  might  not  be  necessary  to  obtain  a  good 
approximation  of  the  optimal  solution.  The  effect  of  the  smooth  components  can  be  captured 
by  a  coarse  model  containing  a  relatively  small  number  of  design  variables  and  thus  can  be 
solved  by  the  MDO  approach  with  a  relatively  low  computational  cost.  This  will  require 
simple  models  for  the  flow  (panel  method  or  small  disturbance  full  potential  on  a  coarse 
grid)  coupled  with  a  plate  model,  or  coarse  finite-element  model,  for  the  structure. 

If  indeed  the  decoupling  condition  holds,  as  discussed,  we  propose  that  the  problem  be 
solved  in  two  stages  as  illustrated  in  Fig.l.  In  the  first  stage,  the  MDO  approach  will  be 
applied  on  a  coarse  model.  The  second  stage  starts  with  the  solution  of  the  MDO  algorithm 
and  the  refined  problem  is  solved  with  the  disciplinary  algorithm,  thus  avoiding  the  enormous 
complexity  of  the  MDO  algorithm  when  applied  on  the  fine  scale  model.  We  claim  that  the 
resulting  design  will  be  a  good  approximation  of  the  optimal  solution.  We  emphasize  that 
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this  is  possible  due  to  the  loose  coupling  between  the  two  disciplines,  otherwise  the  proposed 
approach  will  result  in  a  poor  approximation  of  the  optimal  design  since  the  result  will  retain 
high  residuals  of  the  multidisciplinary  optimality  conditions.  In  that  case  the  MDO  approach 
should  be  applied  also  on  fine  scales  and  a  preconditioner  should  be  used,  as  illustrated  in 
Eq.(5.5),  to  achieve  effective  convergence. 

Finally,  how  far  can  we  extrapolate  the  conclusions  from  this  model  problem  to  a  more 
realistic  model? 

As  for  the  aerodynamic  model,  it  is  shown  in  [12]  that  an  identical  symbol  for  the  aerody¬ 
namic  part  of  the  Hessian  is  obtained  when  using  Euler  equations  instead  of  the  full  potential. 
The  analysis  for  the  Navier-Stokes  equations  has  not  been  completed  yet.  Shocks  were  also 
neglected  in  the  aerodynamic  model,  but  we  postulate  that  they  are  not  going  to  change  the 
main  conclusion  since  shocks  have  a  global  effect  and  are  not  likely  to  affect  the  conditioning 
of  the  Hessian. 

As  for  the  specific  modeling  which  we  have  chosen  to  analyze,  since  there  are  many 
different  models  for  the  cost  function  and  different  constraints  depending  on  the  application, 
it  is  impractical  to  analyze  them  all.  In  the  model  discussed  in  this  paper  we  used  a  penalty 
term  in  the  cost  function  to  account  for  constraints  on  the  structure  deformation.  However, 
if  using  inequality  constraints  instead  of  penalty  terms  it  is  not  clear  how  the  coupling  of  the 
two  disciplines  will  be  affected.  In  practice  most  of  the  constraints  are  not  binding  at  the 
solution,  and  therefore  are  effectively  of  small  number.  When  introduced  in  small  numbers, 
we  anticipate  that  they  are  not  going  to  change  the  basic  structure  of  the  Hessian  near  the 
minimum. 
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COARSE  MODEL  ‘  H 


Figure  1:  An  optimization  strategy  to  solve  aeroelastic  optimization  problems  in  case  of 
decoupling  as  defined  in  Sec. 5.  Apply  the  MDO  approach  on  a  coarse  model  followed  by  a 
disciplinary  serial  approach  on  fine  scales.  The  result  should  be  a  good  approximation  of  the 
multidisciplinary  optimal  solution. 
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